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HOMOGENISATION OF A ROW OF DISLOCATION DIPOLES 
FROM DISCRETE DISLOCATION DYNAMICS* 

STEPHEN JONATHAN CHAPMANt, YANG XIANG* , AND YICHAO ZHU* 


Abstract. Conventional discrete-to-continuum approaches have seen their limitation in de¬ 
scribing the collective behaviour of the multi-polar configurations of dislocations, which are widely 
observed in crystalline materials. The reason is that dislocation dipoles, which play an important 
role in determining the mechanical properties of crystals, often get smeared out when traditional 
homogenisation methods are applied. To address such difficulties, the collective behaviour of a row 
of dislocation dipoles is studied by using matched asymptotic techniques. The discrete-to-continuum 
transition is facilitated by introducing two field variables respectively describing the dislocation pair 
density potential and the dislocation pair width. It is found that the dislocation pair width evolves 
much faster than the pair density. Such hierarchy in evolution time scales enables us to describe 
the dislocation dynamics at the coarse-grained level by an evolution equation for the slowly varying 
variable (the pair density) coupled with an equilibrium equation for the fast varying variable (the 
pair width). The time-scale separation method adopted here paves a way for properly incorporating 
dipole-like (zero net Burgers vector but non-vanishing) dislocation structures, known as the statis¬ 
tically stored dislocations (SSDs) into macroscopic models of crystal plasticity in three dimensions. 
Moreover, the natural transition between different equilibrium patterns found here may also shed 
light on understanding the emergence of the persistent slip bands (PSBs) in fatigue metals induced 
by cyclic loads. 
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1. Introduction. It is well known that the plastic deformation of crystalline 
materials is carried by a large number of atomistic line defects, i.e. dislocations. 
Hence macroscopic models of crystal plasticity can be established by formulating the 
dynamics of many dislocations. As an idealised (but also practically useful) case, the 
dynamics of straight and mutually-parallel dislocations have been intensively studied. 
These translationally invariant dislocations can be treated as “poles” on one of the 
planes perpendicular to all dislocation lines. These straight dislocations, like electrical 
charges, have signs depending on their line directions with respect to the slip direction, 
known as the Burgers vector. Abundant experimental evidence suggests that a good 
understanding of the collective behaviour of many straight dislocations is important 
for controlling the mechanical properties of crystals. One example is found inside a 
single-crystalline fatigued copper specimen induced by cyclic loads [14]. Before the 
saturation point is reached, the inner configuration of the copper specimen takes a 
“channel-vein” structure as shown in Fig. 1(a). A vein consists of many almost straight 
and closely spaced edge dislocations and the veins are separated by channels where 
the dislocation density is relatively low. Beyond the saturation point, a characteristic 
ladder-shape structure (known as persistent slip bands (PSBs)) forms as shown in 
Fig. 1(b). The walls of the ladders also consist of straight edge dislocations. The 
mechanism governing the transition from the channel-vein to PSB structures is still 
unclear, and a study of the collective behaviour of edge dipoles can be of great help 
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(b) PSB structure 


Fig. 1. Dislocation patterns in the early stage of metal fatigue induced by cyclic loads. 


to its understanding. 

One way to reveal the role played by these straight dislocations during the forma¬ 
tion of PSBs, is by using two-dimensional (2D) discrete dislocation dynamical (DDD) 
models, where all dislocations are tracked individually (e.g. [3]). Nevertheless, it is 
still difficult to get a clear idea of the mechanism that governs dislocation pattern for¬ 
mation in crystals from DDD simulations. Hence there is still a need to investigate the 
dynamics of dislocations at the continuum level, where materials substructures are de¬ 
scribed by a dislocation density distribution. In principle, a dislocation-density-based 
continuum model should be obtained through a rigorous averaging of its underlying 
2D DDD model. However, existing discrete-to-continuum approaches struggle to up¬ 
scale multi-polar configurations of straight dislocations. The reason is as follows. At 
room temperature, dislocations (of edge type) are in general constrained in their own 
slip planes. As a result, a positively defined edge dislocation and a nearby negatively 
defined edge dislocation which is not on the same slip plane tend to lock each other by 
forming a pair of dipole rather than to annihilate each other. Since the locking stress 
between the two components of a dislocation dipole scales with the intra-dipolar spac¬ 
ing r by 1/r, a relatively large externally applied stress is needed in order to mobilise 
the constituent dislocations of a dipole. Hence the presence of dislocation dipoles 
may effectively increase the strength of a crystal. When traditional homogenisation 
methods are applied, however, dipoles, despite being crucial in determining the ma¬ 
terial mechanical properties, average to zero and they play no role in the continuum 
approximation. Owing to this, traditional homogenisation techniques are only appli¬ 
cable when investigating the collective behaviour of many dislocations of the same 
sign, i.e. the geometrically necessary dislocations (GNDs) (e.g. [7, 16, 18]). The 
collective behaviour of an arbitrary multipolar configuration of dislocations is only 
considered in a phenomenological or statistical manner [4, 8, 9]. There have also been 
works where each dipole pair is treated as one object so that the traditional homogeni¬ 
sation method can be applied [10]. As shown by our analysis, this only works for the 
case where the slip plane spacing is much smaller than the inter-dipolar spacing. 

Capturing dipole-like structures at the continuum level is also a bottom neck prob¬ 
lem for establishing a three-dimensional (3D) dislocation-density-based continuum 
theory of plasticity. The density distribution of GNDs in 3D space, where disloca¬ 
tions can be curved, is represented by the Nye’s dislocation density tensor [15], which 
only accommodates the gradient of (macroscopic) plastic strains. One missing part is 
the role played by statistically stored dislocations (SSDs), whose physical dimensions 
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are too small to be distinguished from a dislocation-free state in the Nye’s dislocation 
density tensor. Similar as dislocation dipoles discussed above, some SSD structures 
also play a role in determining the (macroscopic) plastic properties of crystals. During 
the past two decades, many valuable works have been done in order to improve the 
framework based on the Nye’s dislocation density tensor (e.g. [1, 2, 5, 6, 11, 13, 17]), 
but the formulation of SSDs at the continuum level is either phenomenological or 
statistical up to date. There are other dislocation configurations that are not prop¬ 
erly included in the framework based on the Nye’s dislocation density tensor, such 
as dislocation interactions with other types of crystalline defects (e.g. Frank-Read 
sources, grain boundaries). Therefore, a pivotal question to be answered for estab¬ 
lishing a solid dislocation-density-based theory of plasticity is, “how should SSDs as 
well as other structures missing in the framework for GNDs be properly formulated 
on a coarse-grained scale?” Part of this question has been answered through the 
establishment of a continuum model of plasticity, where a set of dislocation density 
potential functions (DDPFs) are employed to represent the dislocation substructures 
on a single slip plane [19, 20, 24] and in three-dimensional space [25]. The micro-scale 
mechanisms that are well incorporated into the continuum model characterised by 
DDPFs are the dislocation line tension effect [19], the grain boundary structures [21] 
and the operation of dislocation sources of the Frank-Read type [24]. The hints of 
how to rigorously incorporate SSDs at the continuum level can be found from the 
analysis presented in the current paper. 

Motivated by these issues, the collective behaviour of a row of dislocation dipoles 
is studied here. The discrete-to-continuum transition is facilitated by the introduction 
of two field variables respectively describing the dislocation pair density potential and 
the dislocation pair width. By using asymptotic analysis, we derive coupled evolu¬ 
tion equations for these two field variables. Actually we show that the time scales 
associated with the evolution of the two field variables are different. The dislocation 
pair width, which moves in response to the resolved shear stress at leading order, 
varies on a time scale much shorter than that associated with the dislocation pair 
density, which moves in response to the “stress gradient” (coming from the resolved 
shear stress at the next order). Hence if viewed on the slower scale, fast-varying 
mechanisms take place so quickly that only their steady (or equilibrium) states need 
to be taken into account. As a result, the collective behaviour of a row of disloca¬ 
tion dipoles at the continuum level can be described by an evolution equation for the 
slowly varying dislocation pair density coupled with an equilibrium equation for the 
dislocation pair width. Such discrete-to-continuum approaches asymptotically sepa¬ 
rating active processes according to their associated time scales may pave a way for 
the incorporation of SSDs at the continuum level. Moreover, a transition between 
two distinct dipolar patterns due to instability, which was originally discovered in 
periodically distributed dipoles [23], is also seen here, and the transition may have 
some role to play in understanding the formation of PSBs. 

The paper is arranged as follows. The governing equations for the DDD model, 
which we take as our reference model, are written down in §2. After the introduction 
of the variables needed for the discrete-to-continuum transition in §3, we derive for the 
asymptotic expressions of the resolved stress field in §4. Then the governing equations 
for equilibrium states and the dynamics at the continuum level are presented in §5. In 
§6, the equilibrium states at the continuum level are further analysed and a natural 
transition between different equilibrium patterns is found. In §7, the accuracy and 
efficiency of the derived continuum model are studied. The article concludes with 
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further discussion in §8. 



L 


Fig. 2. The x-y plane is one of the planes perpendicular to all dislocation lines ,and b is 
the Burgers vector. The numbers of positive and negative dislocations are identically N + 1. All 
positively oriented dislocations are located on one slip plane, which degenerates to the x-axis here, 
and all negatively oriented dislocations are put on another slip plane at a distance of s from the x-axis 
(given by y = s). The n-th dislocation pair consists of the n-th positive and negative dislocations, 
whose locations are set to be at (p„, 0) and ( q n ,s ), respectively. The x-coordinate for the center of 
the n-th pair is denoted by x n given by Eg. (3.1). Here the length of the domain of interest L equals 
1 after non-dimensionalisation. 


2 . Dynamics at the level of discrete dislocations. Here we consider the 
case of a single slip system associated with the Burgers vector denoted by b, and all 
dislocations here are straight, mutually parallel and of edge type. The problem is thus 
reduced to one of the planes that are orthogonal to all dislocation tangents. Here the 
plane of interest is set to be the x-y plane as shown in Fig. 2. If we choose b = (b, 0) 
with b > 0, each dislocation can thus be treated as a signed point in x-y plane. Here 
we set a dislocation with its line direction pointing outward the paper plane (see [12] 
for details) to be a “positive dislocation” and denoted by “_L”. A dislocation with its 
line direction pointing inward the paper plane is set to be a “negative dislocation” 
and denoted by “T”. 

The configuration we consider is shown in Fig. 2. There are N + 1 positive 
dislocations lying on the slip plane characterised by the 2 -axis, while TV + 1 negative 
dislocations are put on another slip plane at a distance of s from the 2 -axis. The n-th 
dislocation pair is set to consist of the n-th positive and negative dislocations, which 
are located at (p n , 0) and ( q n ,s ), respectively. 

Concerning dislocation motion, we employ a dislocation mobility law, which only 
allows dislocations (of edge types) to glide within their slip plane at a speed propor¬ 
tional to their on-site resolved shear stress. Under this rule, the motion of the n-th 
positive dislocation is governed by 

(2-1) v+ = ^jj- = m g b(T int (p n , 0) + T e xt (.Pn ; 0)), 

where denotes the speed of the n-th positive dislocation along 2 -direction; Ti nt ( 2 , y ) 
is the internal resolved shear stress field at (x, y) resulting from the dislocation- 
dislocation interactions; r ex t(x, y) denotes the externally applied resolved shear stress 
at (x, y); m g is the dislocation glide coefficient; b = |b|. 

To facilitate further analysis, we consider the problem in a non-dimensional sense, 
that is, all spatial variables are scaled with L; all stress components are scaled with 
fiNb/(2n(l — v)L) and time t is scaled with 2w(l — v)L 2 /(ftm g Nb 2 ), where /i and v 
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are the shear modulus and Poisson’s ratio, respectively. Hence the non-dimensional 
version of Eq. (2.1) becomes 

(2-2) = Ti nt (p n , 0) +T ex t(Pn,0) = T tot (p n ,0), 

at 

where rt 0 t denotes the (non-dimensional) total resolved shear stress field. Similarly, 
the (non-dimensional) gliding speed of the n-th negative dislocation is governed by 


(2.3) 


d q n 
d t. 


Pint (Sint s) 7~ext(Qn> s) — Ptot (*7nj ^) • 


A comparison between Eqs. (2.2) and (2.3) suggests that a positive and a negative 
dislocation move in opposite directions under the same resolved shear stress field. 

The (non-dimensional) internal resolved shear stress field Tj nt is calculated by the 
superposition of the resolved shear stresses due to all individual dislocations [12]: 


(2.4) 

and 

(2.5) 


1 N i | N 

Pint (Pn j 0) = jt J2 __ _ _ . ~ 


1 Qj){iPn - Qj) ~ s ) 


Pint 


(qn,s) = -J2 


N J^o Pn ~ p i N j^o (( Pn ~ q i? + s2 ) 5 


1 (in ~ Pj){{q n - Pj) 2 ~ s 2 ) 1 1 


--Y- 


N fr'o ((?n - Vi? + S 2 ? N p'oQn- q? 


The dynamics at the level of discrete dislocations is thus given by Eqs. (2.2) - (2.5), 
which form a closed system of ordinary differential equations for the 2(iV-|-l) unknowns 
{p4n=0 and {<?n}n=0- 

3. Preparation for discrete-to-continuum transition. Usually the number 
of dislocations in crystals is very large. Hence it is sensible to consider the collective 
behaviour of the system governed by Eqs. (2.2) - (2.5). Mathematically, this can be 
achieved by examining the asymptotic behaviour of the system as N —> oo. The 
expected outcomes are the evolution equations of some continuously defined variables 
that characterise the dislocation substructures. In this section, we will introduce the 
held variables needed for the discrete-to-continuum transition. 

Given a large N, the length scale associated with the discrete dislocation dy¬ 
namical model given by Eqs. (2.2) to (2.5) is the spacing of neighbouring discrete 
dislocations, i.e. 0(1/N), so that individual dislocations can be observed. We now 
want to describe the same dynamical relation by a model associated with an 0(1) 
length scale, where a continuous dislocation density distribution is considered rather 
than isolated dislocations. 

To facilitate such a transition, we first define x n to be the x coordinate for the 
center of the n-th dislocation pair as shown in Fig 2 


(3.1) 




Pn + q n 
2 


Then we introduce a continuous function of (non-dimensional) time and space denoted 
by ((t,x), such that 


C(t, x n ) 


(3.2) 


N 


= Qn ~Pn- 
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Here £ (t,x) is a field variable defined for x £ [0,1], and its value at x n measures 
the width of n-th dislocation pair scaled by N. Since the spacings of neighbouring 
dislocations are 0(1/N), £(t,x) ~ 0(1) as N —> oo. At the continuum level, £ is 
employed to characterise the local pattern of dislocation dipoles. 

Throughout the paper, a subscript n or j affiliated with a field variable such as £ 
indicates that the field is evaluated at x = x n or at x = Xj . respectively; for example, 
£n = £ (t,x n ). Here we consider the case where £ £ [0,1/2], and the properties for 
£ £ (—1/2,0) can be studied likewise. 

From Eqs. (3.1) and (3.2), p n and q n can be expressed in terms of x n and £ n 
respectively by 


(3.3) 


„ _ Cn £n 

Pn — x n 2 jy ’ 9ra — x n + ^ jy ’ 


We now introduce another field variable, the dislocation pair density potential 
<f>(t, x), such that 

Ti 

(3-4) <j> n = <j>(t,x n ) = —; 

this definition is by direct analogy with the dislocation density potential functions 
defined in [19] or [25]. It can be shown by following the same argument presented 
in [19] that the density distribution of the dislocation pairs denoted by p can be 
calculated by p = d<fi/dx. Throughout this paper, the inputs (t,x) for <j> and £ are 
omitted if no ambiguities are caused. Moveover, a dash is added to a variable to 
denote its derivative with respect to x, for example p = <j)'. 

Here we only consider the case when s, the (non-dimensional) spacing between 
the two slip planes, is 0(1/N), as N —> oo. This implies that s can be rescaled by 


(3.5) 


s = 


S_ 

N’ 


where S ~ 0(1). When s ~ 0(1), the interaction between dislocations from different 
slip planes become long-range, and the configurations can be studied by applying 
conventional homogenisation approaches. 

At the continuum level, the dislocation substructures are expected to be described 
by the two field variables <f> and £ and the goal now is to look for their governing 
equations by taking the asymptotic limit N —> oo of Eqs. (2.2) - (2.5). 

4. Asymptotic behaviour of the resolved shear stress field. To accom¬ 
plish the discrete-to-continuum transition we use the following procedure. Given a 
quantity defined in a discrete sense, we first asymptotically express the values at 
(p n , 0) and (q n , s) by functions of x n , for any integer n £ [0, N], In this way the equa¬ 
tions at the discrete level can be transformed into asymptotic equations for <j> and £, 
which only hold at every x n . Then by using the fact that x n is densely distributed 
throughout the whole domain, we replace x n by x to turn the obtained equations to 
corresponding integro-differential equations of <fi and £, which hold for all x. 

Following this strategy, we start by considering the asymptotic behaviour of the 
internal resolved shear stress T m t(p n ,0) and T m t(q n ,s), given by Eq. (2.4) and (2.5), 
respectively as N —> oo. First, an interval is introduced to the n-th dislocation 
pair, such that the x-coordinates of the centers of its 2 K neighbouring pairs fall inside 
as shown in Fig. 3. The number K here satisfies 1 « K <C N. Throughout this 
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Fig. 3. Given the n-th dislocation pair, the x-coordinates of the centers of its 2 K neighbouring 
pairs fall inside an interval, defined to be the inner region whose size is 0(K/N). Mathemat¬ 
ically, this interval is given by Eq. (4.1). The outer region is defined to be the interval into which 
the x-coordinates of the centers of all other dislocation pairs fall. Mathematically, it is given by 
Eq. (4.2). 


paper, defined in this way is termed as the “inner region”. Mathematically, Cl™ n 
is represented by 


(4.1) 


n ■ 




— <*}. 

TV - TV J 


It can be seen that the length of is 0(K/N). Similarly we define the “outer 
region” by 


(4.2) 


o" = 

^ ^out 




TV 


> 


K 

TV 


which contains the x-coordinate of the centers of all other dislocation pairs. Then we 
estimate Pi n t(Pn,0) in Eq. (2.4) by decomposing it into two parts 


(4.3) 


Pint (Pni 0) = r£ t (p n , 0) + T^(p n , 0), 


where T^ t (p n , 0) denotes the resolved shear stress due to all dislocations associated 
with the inner region 


n+K .. 

(4.4) C,(p„,0) = - £ -—-i, £ 


1 (Pn - 1j)((Pn - q 0 ) 2 - S 2 ) 


TV . ^ Vn-Pj TV ^ ((p n - qj) 2 + s 2 ) 2 

j=n—K J j=n—K J 

3+n 


and r/nt* (Pn ; 0) denotes the resolved shear stress due to all dislocations associated with 
O" • 

^ ^OUt ’ 


(4.5) 


out 
' int 


(Pn,0) 


1 

TV 


E 

0<j<n—K 

n-\-K<j<N 


\Pn~ Pj 


(Pn -qj)((Pn ~ qj ) 2 ~ S 2 A 
((Pu-q^+s 2 ) 2 )' 


It is worth noting that the decomposition suggested by Eq. (4.3) only holds for dislo¬ 
cation pairs that are not too close to the boundaries, i.e. K < n < TV — K. 

We will perform an inner and an outer region approximation to calculate the 
asymptotic limit of r™ t (p n ,0) and r-^. t (p n ,0) respectively, as TV —► oo. Then we put 
the results together to get an approximation to Tint (Pn, 0). The asymptotic behaviour 
of Ti nt (q n ,0) as TV — > oo will be studied likewise. 
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4.1. Inner region approximation. In order to get an asymptotic expression 
for T;“ t (p n , 0) as TV —> oo, we first look for the expansion of each term in the summa¬ 
tion appearing in Eq. (4.4). Given a dislocation associated with the inner region, the 
distance from its x-coordinate (for example, pj or qj ), to x n , the x-coordinate of the 
center of n-th dislocation pair is small compared to the length of the computational 
domain (which equals to 1 after non-dimensionalisation). Hence we use Taylor ex¬ 
pansion to asymptotically express pj and qj in terms of x n . This is accomplished in 
two steps: first we relate pj and to Xj and then we relate Xj to x n . The first step 
has been achieved by Eq. (3.3). For the second step, we re-write Eq. (3.4) by 

,, . , j n 7 — n ., . j — n 

( 4 - 6 ) (j){t,x j ) = — = — + J —^- = <j)(t,x n )- 

Applying q i _1 to both sides of Eq. (4.6), noting that < f>~ 1 (j/N) = Xj, gives 

(4.7) xj = (j)- 1 ^(i, x n ) + ■ 


Since < | <1 for all Xj £ we expand Eq. (4.7) in terms of pp to obtain 
(4.8) 


Xj ~ Xji ~\~ 


j-n J_ _ (j - n) 2 


TV 2 


2 (&) 


" , O'- n) 3 (3«) 2 -<C) | c 


TV 3 


6 (^>n ) 5 


A 4 
TV 4 


Recall that an index n on (j) or £ denotes that the evaluation is made at x n . Using 
Eq. (4.8), we Taylor expand Q about x n to give 

/a n\ /■ ru ^ a , 1 U~ n )C , 1 (j'-«) 2 (CnY n -^Cn) , ^ ( K?> 

(«) G = + + w --+ °(vs 

Combining Eqs. (3.3), (4.8) and (4.9), we asymptotically express pj near x n by 
(4.10) 


Pj ~ X n + — • 


j~n Cn 


1 ({j- n)Cn , (3 - n) 2 


27V 2 


(, 3-n ) 2 ( , ti-n)(€) 2 


N 3 

Similarly we find 

(4.11) 


C" 

Sn 


«4 ) 3 

U - n Wn 


qj ~ X n + 


TV 

(.3 - n) 2 
N 3 


4 (^) 3 


j-n Cn 


2 WnY 4(^)2 6 (<#J 4 

(j~n)C (. j~n ) 2 </>" 


+ 0 — 


TV 4 

TV 4 


27V 2 V < W4) 3 

_ (j - n)(j\ 

4(^)3 2(^)5 4 (^) 2 6(^) 4 


c>" u-n)m 


+ o 


/ TV 4 ' 


Then incorporating Eqs. (4.10) and (4.11) into (4.4), we obtain the expansion of 
T k?t(Pn, 0 ) as 


7int(Pn,0) ~ (?r^) ■ G 0 (2n£ n <j >' n , 2ttS^) + 


2 Cn<Pn Cn4> n 


K 


K 2 


~ WIT • Gu(27re„<(4, 2 tt5^) - ^4^ ' 2nScf>' n ) 

N<P n 


TV 




TV 


G 13 (2ir( n cl)' n ,2TrS(l)' n ) + o , 


(4.12) 
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where 

(4.13) 


G 0 (a,p) 


sin a /3 sin a sinh /3 

cosh f3 — cos a (cos a — cosh /3) 2 ’ 


(4.14) 

. 1 a sin a + 2/3 sinh /3 

Gn(a,P) = --- 2(cosa _ cosh/3) 


+ 


+ 


/3 3 sinh/3(l — cos a cosh (3 + sin 2 a) 
4(cos a — cosh /3) 3 


5/3 2 (1 — cos a cosh /3) 3 a/3 sin a sinh f3 

4(cosa — cosh/3) 2 2(cosa — cosh/?) 2 
a/3 2 sin a(l — cos a cosh /3 — sinh 2 (3) 
2 (cos a — cosh / 3) 3 


(4.15) 


Gi2(o;, (3) 


7 ra(l — cos a cosh (3) 
2 (cos a — cosh /3) 2 


7 ra/3 sinh (3(1 — cos a cosh (3 + sin 2 a) 
2(cosa — cosh/3) 3 


and 

(4.16) 

G\s{a,f3) = 

7 rsina / 1 3/3 sinh/3 /3 2 (1 — cos a cosh/3 — sinh 2 f3) 

2 \ cos a —cosh/3 (cos a — cosh/3) 2 (cos a — cosh /3) 3 


The detailed derivation of Eq. (4.12) is given in the supplementary materials. We will 
find that the internal resolved shear stress components accounting for the pair density 
evolution arise at 0(1/N). Thus, unless specified, the expansions to all resolved shear 
stresses will be truncated at o(l/N). To ensure this accuracy, we further choose 
K ~ VN. 

4.2. Outer region approximation. For Xj belonging to the outer region, the 
expansion given by Eq. (4.8) no longer holds since (j — n)/N can grow as large as 
0(1). However, according to Eq. (4.2), we have 


(4.17) 


K 

N 


< 1 4>{ x n) - 


4>'{c 0 )\Xn - Xj |, 


where cq takes some value between xj and x n . Eq. (4.17) suggests that 


(4.18) 


|Xj - x n \ » l/N, 


for all Xj £ f2" ut . Eq. (4.18) implies that the distance between x n and the centre 
of a dislocation pair associated with H" ut is much larger than the spacing between 
neighbouring dislocation pairs. We obtain the expansion of the resolved shear stress 
at (p n , 0) due to the j-th pair in two steps: first, we relate pj and qj to Xj, and relate 
p n and q n to x n by using Eq. (3.3); then we expand as N —> oo using Eq. (4.18). 

Following these two steps, we asymptotically expand the resolved shear stress at 
(pn,0) due to the j-th pair, i.e. the j -th term in the summation of Eq. (4.5), by 

± . ( 1 _ (Pn ~ qj){{Pn - qj ) 2 - ( S/N) 2 ) \ _ l_ (j 

N \Pn - Pj ((Pn - qj ) 2 + ( S/N ) 2 ) 2 J N 2 [x n - Xj ) 2 

1 CnQ-3^ 2 1 3C J -(65 2 -Cg) + 18S 2 C»-Cj 

N 3 ( x n -Xj ) 3 N 4 A{xj-x n y \K 3 )' 


(4.19) 
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We see from Eq. (4.19) that the leading-order effect of the stress at (p n , 0) due to 
the positive dislocation at (pj, 0) cancels with that due to its pair partner located at 

It is worth noting that, when the truncation is made, 1 /|Xj — x n \ can be as 
large as 0(N/K). Besides, the summation made for outer region approximation 
involves almost N terms. Therefore, to ensure an accuracy of o(l/N) for a resolved 
stress component, the truncation made at each term in the summation should be at 
o(l/N 2 ). This is the reason that a truncation at 0(1/K 3 ) is made in Eq. (4.19), given 
K ~ VN. 

Incorporating Eq. (4.19) into (4.5) gives the expansion of (p n , 0): 


(4.20) 


out 

T int 


(Pn,0) 


1 

JP 


E 

0<j<n—K 

n-\-K<j<N 



l 

W 


E 

0 <j<n—K 
n-\-K<j<.N 


CnCj ~ 3 S 2 

(X n - Xj ) 3 


J_. V 3 Cj(65 2 -Cg) + 18g 2 Cn-Cj | n f N \ 
TV 4 4(xj — x n ) 4 \K 5 J 

o <j<n-K v 3 n ' y / 

n-\-K<j<N 


To evaluate the summations appearing in Eq. (4.20), we make use of the Euler- 
Maclaurin formula. The details are in the supplementary materials, and the result is 


(4.21) 


out 
' int 


(Pn, 0) 


2 1 i / 

’ (Cn0n) 4" 2 * C n4*n "b ’ ( 

1 r*" (^(q)C(a))'da / 1 \ 

NJx o X n -a \N J 


#»Co 

Xn ^0 


where the integral is evaluated in the sense of principal value. 


_^nCn_\ 
Xn Xn J 


4.3. Total resolved shear stress. Now we put the results from the inner ex¬ 
pansion by Eq. (4.12) and from the outer expansion by (4.21) together to obtain the 
expansion of r in t(p n ,0) as 
(4.22) 


Tint (Pm 0) 


(t</ 4) • Go(27T</4C„, 2 TrS(/)' n ) + ^ • 


$.Co 


_^nCn_\ 

Xn Xn J 


| 1 p W(a)aa))'da < 
NJ X0 x n -a N<t>' n 

MEZ. Gl2{27r ^ c „, 2 ^) 


Gii(2n<l>nCn,2nS<j>' n ) 

- . Gl3 (27r^Cn, 2 tt 5^) + o (E) . 


It is worth noting that the 0(1/K) and 0(1/K 2 ) terms from the inner expansion 
cancel with their counterparts from the outer expansion. As a result, no trace of the 
intermediate parameter K is seen in Eq. (4.22). 

The external stress r ex t at (p n , 0) can also be expanded near (x n ,0) 

( 4 -23) Text [jPn ? 0) ~ T° xt (x n ) - E ^ “ + O ^ > 


where for ease of notation we have written T® xt (x) = r ext (x, 0) and dr® xt (x)/dx = 
dr ex t{x , 0)/ dx. Similarly, 


(4.24) T ex t(i?n,-5) 



3T e °xt(zn) , 0 5r e ° xt ( x n )\ , „( 1 ^ 
dx +b ' dy ) +U {n1)’ 
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where dr^ xt {x)/dy = 9r ext (x, 0 )/dy. Here dr® xt /dx and dr® xt /dy are the stress gra¬ 
dients, which capture the difference in the externally applied stress field evaluated at 
each component of a dipole pair. 

Thus the (non-dimensional) total resolved shear stress at (p„,0) is given by 


Ttot (Pn, 0) ~ (t T<j>'n) ■ G 0 (2-JT(j) , n C„, 2? xS<j)' n ) + T e ° xt (l„) 


+ 


N 


<(>oCo 


4 >'n(n 


Xn X 0 Xyi %N 


) + jiL 


(</>'(a)C (o))'da 


(4.25) 




Gu(2n<j> n C n , 2nS<f)' n ) - ‘ G 12 (2n<p' n {n, 2tt5^) 


' c 

nSn 


N 


Gi3(2ir4l n (; n ,2itS4> n ) - 


N 

C n dr° xt (x n ) 
2 N dx 


o — 
1 N 


where we recall that Go, Gu, G 12 and G 13 are defined by Eqs. (4.13) - (4.16). Simi¬ 
larly, 

(4.26) 

rtot{q n , s ) ~ (tt<) • G 0 (2n<t>' n ( n , 2nSfiJ + r e ° xt (x n ) 

+ 1 . ( _ PnSn \ + 1 f XN ^(a)ga)yda 

N \ X n Xq X n X}\j J Nj xQ X n CL 


G n (27r<Cn, 27 rScfr'J + • G 12 ( 2 ^<C n , 2tt S<f>' n ) 


N<f)' n v rn ' N 

finCn /r} jJ /■ \ 1 Cn ^^ext ip^n) 

— ■ G 13 ( 2^„C„, 27t5^„) + —-^- 


5* <9r e 0 xt (:r n ) 


TV 





5. Dislocation dynamical model at the continuum level. In this section, 
we derive the governing equations for the two field variables <j> and £ at the continuum 
level. We first consider the continuous description for the equilibrium states, where 
the total resolved shear stress at each dislocation vanishes. 

5.1. Governing equations for the equilibrium states. When all dipoles are 
in equilibrium, the total resolved shear stress T to t(Pn, 0) and T to t(q n , s) should be zero 
for all n according to the laws of motion (2.2) and (2.3). It is worth noting that since 
the resulting equations are established in an asymptotic sense, we also need to expand 
(j) and ( as 

(5.1) 4>(t,x) ~ g°\t,x) + -—H-• 

and 

(5.2) C(t^)-C ( 0 ) ( i^ ) + ^|l^ + ---, 

respectively. Substituting the above expansions into Eqs. (4.25) and (4.26) gives 


Ttot(Pn,0)~ (t rO^y) • Go(27r(^ ) ) , Ci 0) ,27rS'(4 0) )')+ 7- e 0 xt (a;n) 


J_ ( n _ n dr° xt (x n ) \ 

Ny a b 2 dx J 


+ o 



(5.3) 
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and 


Ttot(q n ,s) ~ • Go(27r(^ 0) ) , C^ ) ,27r5'(4 0) )')+ r e ° xt (x n ; 

^ 5 ' 4 ^ 1 ( d 0) dr° xt (x n ) 0 5r e ° xt ( x n ) 


+ 


N 


dx 


+ S- 


dx 


respectively, where 
(5.5) 


(CrC ttWrcW , r* (wc(a) (o) ydt 


^0 X N 

+c 


+f 

^ Xn 


x n -t 


d.1 dGo(2n((j)n^yCn \ 2nS(<J>n^y) , ^ (1)v SG 0 (27r(<Ai 0) ) , d° ) , 27rS(^ 0) )') 

" -- + ( *>-^ - ' 


T& ” = ’ Gii(27r(^ o) ) / e ) ,27r5(4 o) )0 

\ 4 >n y 

+ ((^yefy • G 12 (27r(4°))'cy 3) ,2^(^°))') 
+ (4 0) )'(Cf)' • G 13 (27r(4°))'cy )) ,2^(4 0 ))'). 


Now letting the right hand side of Eqs. (5.3) and (5.4) vanish and equating coef¬ 
ficients of the same powers of IV, we obtain at leading order, 

(5.7) ^(^y • G , o(27r(</>y )) ) , C^ cl) , 2T:S{(j)^)') + T® xt (x n ) = 0. 

There are N + 1 equations for the 2(7V + 1) unknowns {(f’ndn^o and {Cn°^}^Lo- 

To close the system, we need to proceed to higher order in the expansion. At 
0(1/N), we find 


(5.8) 
and 

(5.9) 


_^n _n 
T a - T b 


d° } dr° xt (x„) 


dx 


= 0 


d° ] dT° xt (x n ) clr e 0 xt (x n ) 


dx 


+ S- 


dx 


= 0. 


Now there are 4(iV+l) unknowns {0n O) }OL o , {Cn 0) }n=o and {Cn^j^Lo and 

3(A r + 1) equations. However, we subtract Eq. (5.8) from (5.9) to eliminate (l 1 * and 
(j)n \ both of which only appear in t", to give 


(5.10) 2r b n + (i 0) ■ + S ■ dT ^ Xn ^> - 0. 

Eqs. (5.7) and (5.10) form a system consisting of 2(iV + 1) equations for 2 (IV + 
1) unknowns ({bn^jOLo and {Cn°' 1 j-^-o)- Henceforth we drop the superscript “(°)”, 
because only the leading-order effects are taken into account. 

Since x n is densely distributed in the domain, we rewrite our equations valid at 
every x n as equations valid for all x. Therefore, we drop the index n and re-write 
Eqs. (5.7) and (5.10) as 


Tt(j)' sin(27T0 , C) / 2t:S cj)' sinh(27r5 , ^) / ) \ 0 

cosh(27r5^>') — cos(27 v<j)'Q \ cosh(27rS<^') — cos(27r<//£) J ext 
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26" 

0 = -Jj- • Gh(27T(/ ) , C, 2 TrS<t>') + 2(0'C)' • Gi 2 (27r<£'C, 2 ttS^) 

(5 - Ub) 9r o 

+ 2(0 / C / ) • G 13 (2^'C, 2irS4>') + 

respectively, where we recall that Gn, G 12 and G 13 are defined by Eqs. (4.14) - 
(4.16). Eqs. (5.11a) and (5.11b) are the two equations for the two field variables </> 
and C, derived at the continuum level when the row of dipoles rest in their equilibrium 
states. It is worth noting that Eq. (5.11a) comes from the leading-order force balance 
and Eq. (5.11b) comes from the difference in the force balance equations obtained at 
the next order. 


5.2. Governing equations for the dynamics. Now we consider reformulating 
the discrete dislocation dynamics governed by Eqs. (2.2) to (2.5) at the continuum 
level by looking for evolution equations for <j> and £. 

We know by definition that </)(t, x n (t )) = n/N at any time t. Taking the derivative 
with respect to t on both sides gives 


(5.12) 


d<j) n d^n 9(j) n 
dt df dx 


According to the definition of x n given by Eq. (3.1), we have 


(5.13) 


d-X'r, _ Vtof (Pn , 0) Dot. {(jn, a) 

~df ~ 2 


where the laws of motion (2.2) and (2.3) are employed. With the asymptotic expan¬ 
sions for Ttot(Pn)O) and T tot (q n ,s) given by Eqs (4.25) and (4.26), respectively, we 
incorporate Eq. (5.13) into (5.12) to get 
(5.14) 


d^n _ 
dt N 
1 

~~ N 


f Gupw^Sn^nSfa) + «Cn) , G 12 (27r^Cn,2^) 


d(j)n 

dx 


«C)Gi3(27r^Cn,27r5<) + 


C^r° xt ( x n ) S dT° xt {x n )\ dcj>. 


dx 


dy 


dx 


Again we drop the subscript n to rewrite Eq. (5.14) as a differential equation valid 
for all x by 


(5.15) 

where 


d<j) 1 


1 

( 

9r e ° xt 

S 

<9r e ° xt \ 

dcj) 

( 1 \ 

~N 1 

( ,l+ 2 ' 

dx 

+ 2 ' 

dy ) 

dx 

-0 U) 


(5.16) n = 2ir S4>')+((f>'(yG 12 (2w(j)X, 2 tt Scf)')+(t>'£'G 13 (2ttc))'£, 2 r xS(j)'). 


Eq. (5.15) can be considered as the evolution equation for (j). 

It can be seen from Eq. (5.15) that the evolution speed of (j) is as small as 0(1/N). 
This suggests that the natural time scale associated with the evolution of <j>, the 
dislocation pair density potential, is characterised by a slow-varying temporal variable 
given by t s = Nt. Eq. (5.15) then gives at leading order 


(5.17) 


fty _ ( , C d^ext , S dr° xt \ 

dt s V 2 dx 2 dy J dx 
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On the other hand, according to the definition of ( by Eq. (3.2), we have 

o\ ^Cn d X n dCn d(^ n (tj x(t)) j\r (d Qn dp n \ 

dt d t dx df \ dt d t ) 


Combining Eqs. (2.2), (2.3), (4.25), (4.26), (5.13) and (5.18) then dropping the sub¬ 
script n, we find 
(5.19) 

d( 2Nir(t>' sin(27n//(C) / 2nS(j)' sinh(27r5 , ^ / ) \ 0 , 

dt cosh(27rS , </>') — cos( 27T(/)'() \cosh(27rS'0') — cos(27T0'^) / ext 

It is seen from Eq. (5.19) that £ evolves as fast as O(N). This means £ should be 
studied at a fast temporal scale characterised by t{ = t/N. Then the leading-order 
equation for £ is 


(5.20) 


27 T(j)' sin(27r<//£) 


ac = __ 

dtf cosh(27r S(f)') — cos(27t^/£) 


1 - 


27t S(f>' sinh(27r5'<^ / ) 
cosh(27r5 , ^ / ) — cos( 27r<^'£) 


— 2 r 

^ / f-i-i 


A comparison between Eqs. (5.17) and (5.20) shows that the evolution of £ is much 
faster than that of </>. Hence £ can be considered varying quasi-static ally on the time 
scale characterised by t s , on which <j> naturally evolves, provided stable equilibria 
exist for Eq. (5.20). In fact, Eq. (5.20) can be written by ^ where T is the 

generalised free energy density with respect to given by 


(5.21) T = log (cosh(27rS'<(/) — cos(27r^')) 


2n(f>' S sinh(27r S<f>') 
cosh(27r S(j)') — cos(27 rC,(j)') 


+ 2 Cr e 0 xt 


Since (j)' is assumed static on the fast scale, the stable equilibria of Eq. (5.20) are 
identified wherever J 7 attains its local minimum with respect to (,. It will be shown 
numerically later that given (f>' and S , stable equilibria exist for Eq. (5.20) when |r,? xt | 
falls below some critical value. 

Therefore, the dynamics of a row of dislocation dipoles at the continuum level 
can be described by the following coupled equations: 

irfi sin(27n//Q / 2nS(j) 1 sinh^nSft) \ 0 

cosh(27r50') — cos(27r0 / C) \ cosh(27r50') — cos(27r </>'£)) ext 


(5.22b) 


^__( T , C && t S dr° xt \ <M =() 

dt s \ b 2 dx 2 dy ) dx 


where n was defined by Eq. (5.16), provided the stable equilibria of Eq. (5.22a) exist. 
Noted that Eq. (5.22a) is effectively the leading order force balance equation (5.11a). 

6. Equilibria at the continuum level. In this section, we will analyse the 
equilibrium states at the continuum level determined by Eqs. (5.11a) and (5.11b). We 
will begin with the case where the externally applied stress vanishes on y = 0. In this 
case, two types of possibly stable configurations are found as a result of the leading- 
order force balance equation and a natural transition between different equilibrium 
patterns due to instability is seen. At the next order, the detailed equations for (f> and 
( corresponding to various equilibrium states will be derived. The analytical results 
will then be validated through comparison with the numerical solutions to the same 
problem by using the DDD model. In the end of this section, we will analyse the 
equilibrium under arbitrary externally-applied stresses. 
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6.1. Equilibria under an external stress field which vanishes on y = 0. 

We now analyse Eq. (5.11a) and (5.11b) by starting with a simple case where the 
externally applied resolved shear stress vanishes on y = 0, i.e. r° xt = 0. Note that 
the stress gradient dr® xt /dy need not vanish. 

6.1.1. Implication from the leading-order force balance equation. When 
r ext = 0; the leading-order force balance equation (5.11a) becomes 

_ sin(27nft'C) _/ 27 t 5 i <// sinh(27r5'(/>') \ 

cosh(27rS'</') — cos(27T0'C) \ cosh(27rS'</') — cos(27 T(p'Q J 

Eq. (6.1) can be regarded as an implicit relation between the two quantities (j)' C, and 
<j>'S. In fact, these two quantities are physically meaningful. Since the pair density <f>' 
can be approximated by the reciprocal of the spacing between two neighbouring pair 
centers scaled by N, and S is the slip plane gap rescaled by N, < j>'S captures the ratio 
of slip plane gap to the inter-spacing of neighbouring pairs. Also since Q/N measures 
the pair width at x according to Eq. (3.2), <j>'£ measures the ratio of the pair width 
to the spacing of neighbouring pairs. 

From Eq. (6.1), there are three possible choices for </ as a function of cj>’ and other 
parameters. 

• Equilibrium Type I when £ = 0. The dislocation substructure is shown 
in Fig. 4(a). Within each dislocation pair, the positive and the negative 
dislocations are vertically aligned. 

• Equilibrium Type II when <(/£ = 1/2. The dislocation substructure is shown 
in Fig. 4(b). Since 4>'C represents the ratio of pair width to pair center spacing, 
c t>'(, = 1/2 suggests that every negative dislocation lies roughly in the middle 
of its two neighbouring positive dislocations. We term the equilibrium of this 
type as a “non-localised” structure, because each dislocation is “shared” by 
its two neighbours. 

• Equilibrium Type III when 

(6.2) ( = \ cos -1 (cosh(2-K Sc/>') — 2-kSi j>' sinh(27r Scj)')). 

27 T(j>' 

The dislocation substructure is shown in Fig 4(c). A positive dislocation here 
is bonded with a negative one to form a real dipole, and the equilibrium of 
this type is named as a “localised structure”. 

It is worth noting that Eq. (6.2) only holds when 

(6.3) —1 < cosh(27rS'<//) — 2nS(j)' srah(2-irS(j)') < 1, 
which numerically gives rise to a range for S(p': 

(6.4) 0 < S(j/ < 0.2465. 

Hence the emergence of Equilibrium Type III is conditional. 

If we set X = </)'(, and Y = (j)'S, the configuration is equivalent to a row of dipoles 
periodic in X 1 which have been studied in [23]. Thus the conclusion regarding the 
stability of the obtained three types of equilibria can be drawn by employing the same 
arguments proposed by [23]: 

• Equilibrium Type I (£ = 0) is always unstable. 

• Equilibrium Type II (<//(/ = 1/2) is only stable when Equilibrium Type III 
does not exist. 
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(a) Type I 


(b) Type II 


(c) Type III 


Fig. 4. Three types of equilibria: (a) £ = 0; (b) /V = 1/2 with non-localised structures formed; 

(c) £ satisfies Eq. (6.2) and localised structures are formed. 


• Equilibrium Type III (£ satisfies Eq. (6.2)) is always stable as long as it exists. 

Another way to investigate the stability of the obtained equilibrium states is to 
look for the local minimum of the free energy density J- with respect to </. When 
Text = 0, T given by Eq. (5.21) are drawn against ( for different <f>'S as shown in 
Fig. 5. It is seen from Fig. 5(a) that when condition (6.4) is not satisfied, there are 



Fig. 5. A stable equilibrium state should correspond to a local minimum of the generalised free 
energy density T given by Eq. (5.21) and r ex t = 0 with respect to f. (a) If <p'S is larger than 0.2465, 
only two types of equilibria exist and Type II is the stable configuration, (b) If 0 < S < 0.2465, a 
transition in stability from Type II to Type III takes place. 


two equilibrium states, and Equilibrium Type II is the stable one. When condition 
(6.4) is met, we have three equilibrium states as shown in Fig. 5, and Equilibrium 
Type III is the stable one. 

Here a natural transition from a non-localised structure (Type II) to a localised 
structure (Type III) takes place as the slip plane spacing gets narrower or equivalently, 
as the pair density decreases. Such a transition may be indicative of the formation of 
the persistent slip bands; further discussion on this issue will be made in § 8 . 2 . 


6.1.2. First-order force balance equation for Equilibrium Type II. Based 
on the solutions to the leading-order equation (5.11a), we now investigate the first- 
order force balance equation (5.11b). Here only stable configurations, i.e. Equilibrium 
Type II and III, are considered. 

When cf>'( = 1/2 (Type II), one can make use of the fact that (//£)' = 0 and 
sin(27r^'C) = 0. This suggests that the terms associated with G\i and G 13 in 
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Eq. (5.11b) both vanish. Therefore, the equation for qb' can be obtained as 

0 = —- • Gn(7T, 2nb<p ) + b —— 
ft dy 

(6.5) = — —- — AirSft' tauhtnSft ) + 5TT 2 S 2 ftft'sech 2 (nSft) 

ft 

— 2n 3 S 3 (ft) 2 ft'sech 2 (irSft) tanh(7r Sft) + S ^° xt . 

Eq. (6.5) is a differential equation for ft, the (non-dimensional) pair density. Its 
solution describes the pair density distribution in equilibrium when all dipoles form 
non-local structures as shown in Fig. 4(b). 

To justify our results for ft and £ calculated from the continuum model, we also 
consider the equilibrium states obtained by the discrete dislocation dynamical model. 
To do that, we simply put N + 1 pairs of dipoles in the domain [0,1] and let the 
system evolve to the steady state following Eqs. (2.2) - (2.5). 

For all the simulation results presented in this paper, we lock one dislocation 
at each end. For example, at the left boundary, we set po = 0. There is no strict 
requirement for qo, except that qo > 0. Similarly at the right end, we let qN = 1 and 
p N < 1. By doing this, the total number of dislocation pairs are conserved during the 
simulation. Correspondingly at the continuum level, this condition is translated by 

(6.6) f ft (t, x)dx = (j){t, 1) — 0) = 1. 

Jo 

The temporal derivatives needed for DDD simulations are approximated by using 
the Euler scheme with time step At<jis chosen by Atdis = 0.025/A. 

Now we compare the results for Equilibrium Type II obtained from the continuum 
and DDD models. For simplicity, we consider the case when dr® xt /dy is a constant. 
Thus we integrate Eq. (6.5) on both sides to obtain 

(6.7) io g ( c ° s i 1 ^ s ) ) + ( C o s ;f4-s) ) +WStal(./S) = (;-^.S>, 

where C is a constant to be determined by condition (6.6). 

We begin with the case when no stress gradient is applied to the system, i.e. 
dr® xt /dy = 0. In this case, Eq. (6.7) suggests ft = 1 and we then obtain ( = 1/2. 
This means that in the absence of applied stress gradients, all dipoles are uniformly 
distributed and the dipoles form non-localised structures suggested by the continuum 
model. To see Equilibrium Type II from the DDD model, one needs S > 0.2465/ ft 
and S is chosen to be 0.3 here. 

Note that in the DDD model, the pair density is approximated by pdis [iPn + 
q n )/ 2) = \/{N(jp n+1 -p n )), and C is approximated by Cdis((p n + 9n)/2) = N{q n -p n ). 

A comparison of the values of ft and C, from the discrete and the continuum models 
is shown in Fig. 6 and good agreement between the two models is seen except near 
the boundaries. There is a boundary layer near each end, where the results from the 
continuum model deviate from its DDD counterpart. This is because the symmetry 
required for the setting up of the inner region 120, given by Eq. (4.1) breaks down. 
However, the goal of this paper is to formulate the collective behaviour of dislocation 
dipoles in the (relatively vast) interior region. It is suggested by the numerical results 
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DDD 

-Continuum 


0.2 0.4 0.6 0.8 

X 

(a) Pair density 



0 0.2 0.4 0.6 0.8 1 

X 

(b) Local pattern 


Fig. 6. Comparison of the pair density and the pair width Equilibrium Type II with the results 
from the discrete dislocation dynamical models in the absence of applied stresses or stress gradients. 
When S = 0.3, dipoles form Equilibrium Type II. Here N = 50. The dipoles take a uniform 
distribution within [0,1]. 


shown below that the influence cast by the boundary layers over the accuracy of the 
continuum approximation in the interior region is limited. Hence the incorporation 
of boundary layers into the continuum framework will be discussed in future work. 

With a non-vanishing stress-gradient, for example, dr® xt /dy = 1, one can again 
calculate (j>' and £ with reference to Eq. (6.7). The results from the two models are 
compared in Fig. 7 and excellent agreement in the interior region is again seen away 
from the two ends. 



(a) Pair density 


DDD 

-Continuum 


0 0.2 0.4 0.6 0.8 1 

X 

(b) Local pattern 


Fig. 7. When the system is applied an stress gradient dr® xt /dy = 1, the dipoles of Equilibrium 
Type II are seen piling-up against the left boundary. Here S = 0.3 and N = 50. 


6.1.3. First-order force balance equation for Equilibrium Type III. Sim¬ 
ilarly, we study the first-order equation (5.11b), when all dipoles are in Equilibrium 
Type III, i.e. condition (6.4) is met. It is recalled that Equilibrium Type III only 
appears for small <p' 5, we consider the asymptotic behaviour of the above equation 
as S —> 0 for simplicity. Thus one can asymptotically solve Eq. (6.2) to get 


( 6 . 8 ) 


S + 


2{irft) 2 S 3 

3 


Eq. (6.8) implies that in this case the pair width is almost the same as the slip 
plane gap. When these two quantities are identical, we call the resulting dislocation 
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structure a 45° dipole. In fact, a 45° dipole is the stable configuration of an isolated 
pair of dipole. We see from Eq. (6.8) that when the two slip planes get close to 
each other (as S —» 0), the mutual interaction between the pair partners becomes 
dominant over the stresses due to all other dislocations, and the dipoles behave as 
isolated dipolar pairs. Incorporating Eq. (6.8) into the first-order equation (5.11b), 
we asymptotically derive an equation for the pair density <fi' in the limit that S —>■ 0 
as 

f)r° 

(6.9) 2-k 2 S 2 cf> 4>" + S ■ = 0. 

dy 

Eqs. (6.8) and (6.9) are valid only when S —> 0. Now we compare their solutions to 
DDD simulation results to show that they can be used as the governing equations for 
many dipoles in equilibrium of Type III at the continuum level. 

Here we still consider the case when dr® xt /dy is constant for simplicity. Hence 
the pair density distribution ((/ can be solved from Eq. (6.9) 

where C is determined by boundary condition (6.6). 

We first investigate the case with no applied stress gradient. From Eq. (6.10), 
we obtain </>’ = 1 and ( is calculated to be 0.1066 from Eq. (6.8). We then compare 
these results with that from the DDD simulations in Fig. 8. Excellent agreement 


DDD 

-Continuum 


0.2 0.4 0.6 0.8 

X 

(a) Pair density 



(b) Local pattern 


Fig. 8. Comparison of results from the continuum and the DDD models for the case where 
there is no applied stress gradient. When S = 0.1, the continuum model suggests that the system 
takes the equilibrium state of Type III with <f> « 1 and f £3 0.1066. Here N = 50. 

between the two models is seen. Here we find again that in the absence of applied 
stress gradient, the dipoles are uniformly distributed. 

Now we consider a non-vanishing applied stress gradient set to be dr° xt /dy = 1. 
By using Eqs. (6.8) and (6.9), we plot </>' and ( against x in Fig. 9 and they are shown 
agreeing well with the outcomes from the underlying DDD model. The comparison 
results shown above suggest that we can use Eqs. (6.8) and (6.9) to describe the 
collective behaviour of a row of dislocation dipoles in equilibrium of Type III. 

6.1.4. Equilibria of mixed types. According to Eq. (6.4), (f>'S = 0.2465 char¬ 
acterises the transition between Equilibrium Type II and III. Therefore, when the 
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(a) Pair density 


DDD 

-Continuum 


0.2 0.4 0.6 0.8 

X 

(b) Local pattern 


Fig. 9. Dipoles of Equilibrium Type III are found piling-up against an applied stress gradient 
to the left boundary. Here S = 0.1, <9r° /dy = 1 and N = 50. 


value of (f)' — 0.2465/S' changes its sign, there should be a change in equilibrium pat¬ 
terns as suggested by the continuum model. This is actually observed in Fig. 10, 
where S is set to be 0.24 and N is chosen to be 100. It is seen from Fig. 10 that the 
dipoles take Equilibrium Type II near the left boundary, and a transition from Type 
II to III is found taking place away from the left end. The continuum model suggests 
that the transition should happen when <// = 0.2465/S « 1.03, which gives rise to the 
dashed line in Fig. 10. It can be checked that Equilibrium Type III roughly emerges 
where <f>' drops below the dashed line. In Fig. 10, it can also be seen that the values 
of the pair density agree well for both equilibrium types, while there is roughly a 10% 
variance in £ for Equilibrium Type II with the change of equilibrium type not so easily 
determined. We will see later that increasing N will bring down the deviation in <f>' 
and </ between the continuum and the DDD models. 




Fig. 10. When S = 0.24, Equilibrium Type II and III are found co-exist. Near the left boundary, 
the dipoles take the equilibrium of Type II. A natural transition from Type II to III is seen roughly 
where the pair density drops below the dashed line characterised by <f>' « 1.03. Here N = 100. 


6.1.5. Summary. To summarise, a row of dipoles may form two types of stable 
equilibria if the applied stress vanishes on y = 0. When <f>'S > 0.2465, the resulting 
equations at the continuum level of the pair density <fi' and (rescaled) pair width £ 
are derived to be Eq. (6.5) and ( = 1/(2</'). When 0 < (j>'S < 0.2465, the collective 
behaviour of a row of dipoles can be approximately described by Eqs. (6.8) and (6.9). 
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6.2. Equilibria under arbitrary externally-applied stresses. Now we gen¬ 
eralise our discussion to the case where the leading order of the external resolved 
shear stress is non-vanishing, i.e. Tg Xt ~ 0( 1). In this case, Eq. (5.11a) may not be 
solved explicitly. However, some analysis can still be done to understand the resulting 
equilibrium configurations. 

If we use the expression for Go defined by Eq. (4.13), we rewrite Eq. (5.11a) as 
(6.11) G 0 (27n//C,27r<//S) + ^ = 0. 

Eq. (6.11) describes the inter-relation of three quantities, ((/)', S<ft and r° xt / 0 ' and we 
define X = ((/)', Y = Sc/)', and T = T° xt /(n^>') to facilitate further analysis. As dis¬ 
cussed in §6.1.1, X and Y measure respectively the pair width and the slip plane gap, 
both scaled by the spacing between the neighbouring dipolar centers. Thus Eq. (6.11) 
can be written as — Go(27tX, 27tF) = T, which suggests that the inter-relation be¬ 
tween X and Y for a given T can be visualised by the contours of —Go(27rX, 27 tF) as 
shown in Fig. 11. It can be observed that on each contour, there exists a Y* (attained 



Fig. 11. Given any Y = T® xt /( 7T0'), a pair of (X. Y), which satisfies Eq. (6.11) should sit on the 
contour —Go{2iuf>'C, ) 2tr(j>'S) with height T. For each Y, there exists a Y* (attained at X* say) such 
that Y < Y*. The locus of such (X *, V* ) lies on the dashed curve. For any Y < Y* (under a given 
Y ), there are two possible values for X. Only those ( X , Y) falling in the shaded region correspond 
to stable configurations. 


at X* say) such that Y < Y*, and the locus of such (X*,Y*) sits on the dashed 
curve in Fig. 11. This means the solution £ to Eq. (6.11) conditionally exists. Given 
Y = S(f >', the solution for X = C,4>' satisfying Eq. (6.11) exists for 

(6.12) Text = 7T0'T < 7T0' • |G 0 (2 ttX*, 2^F*)I. 

The physical interpretation of Eq. (6.12) is that a dipole breaks down to two monopoles 
when the external stress is large. 

It is also observed from Fig. 11 that there exist two choices for X when Y < Y*. 
One way to identify the stability of the candidate solutions is by investigating the 
local minima of the generalised free energy density T by Eq. (5.21) with respect to 
Here we find that the larger one gives rise to a stable equilibrium state after checking 
with the numerical results to be shown later. Hence we conclude that only those 
(X,Y) falling into the shaded region in Fig. 11 correspond to stable configurations. 
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It is worth noting that T > 0 is considered in the analysis presented above. When 
T < 0, we simply let X < 0 and same conclusion will be drawn. 

The above analysis provides us some insight to the equilibrium configurations 
under an arbitrary externally-applied stress. Nevertheless, to find ( and <f>' satisfying 
Eq. (5.11a) and (5.11b), one has to turn to numerical methods. 

7. Comparison of the continuum model with its underlying DDD model. 

Now we compare the simulation results obtained by applying the continuum model 
and the DDD model to same dynamical processes. For simulations at the discrete 
level, the set-up and procedure is as in §6.1.2. To numerically implement the contin¬ 
uum model, we discretise Eqs. (5.22a) and (5.22b) with step Ax in space and At con in 
time. At each time step, we use the following procedure to update the two variables 
(j> and £. With (j> computed from the previous step, we (numerically) solve Eq. (5.22a) 
to update the value for ( at each spatial grid point. It is worth noting that following 
the analysis in §6.2, we need to ensure the computed ( is associated with a stable 
equilibrium state. Then we use Eq. (5.22b) to update <t>. For the simulation results 
presented here, At con was chosen to be 1.25Ax 2 . 

Our goal here is to check the accuracy and the efficiency of the continuum model 
with reference to its underlying DDD model. To measure accuracy, we define 

(7.1) Err,y=max —- 

x£l p dis 

where p d i s denotes the density computed by the DDD simulations; we choose I = 
[0.1,0.9] to avoid the inherent difference between the two methods near the two 
boundaries. Thus Err ^ is used as a measurement of the relative error of the pair 
density caused by the discrete-to-continuum transition. In a similar sense, we define 
a measurement of the relative error of the pair width by 

(7.2) Err<* = max-— 

X ^I Cdis 

The parameters chosen for the first set of numerical examples are S = 0.3, N = 50, 
T ext = 0-5 and dr® xt /dy = 1. In Table. 1, Eny and Err^ at various times are listed. 
Note that the time t in Table 1 is measured in unit 2-7 t(1 — v)L 2 /(m g pb 2 ) with L 


t 

1 

2 

5 

10 

20 

26.4 

Err 0 / 

0.0150 

0.0117 

0.0088 

0.0077 

0.0079 

0.0079 

Err ? 

0.0797 

0.0801 

0.0810 

0.0815 

0.0818 

0.0818 


Table 1 


Defined by Eq. (7.1), Err^ provides a measurement of the relative error of the pair density 
caused by the discrete-to-continuum transition. Similarly Err £ given by Eq. (7.2) provides a mea¬ 
surement of the relative error of the pair width f/7V. Here S = 0.3, r® xt = 0.5, dr® xt /dy = 1 and 
N = 50. Here t is measured in unit 27t(1 — v)L 2 /(m g pb 2 ). Simulations by the two models both stop 
at t = 26.4, when the difference in the dislocation positions in DDD simulations between this and 
the previous time step is no more than 10 — ^At^is- Err^r and Err £ are listed at various times. 

recalled to be the computational domain size. The simulations based on both the 
continuum and DDD models are stopped at t = 26.4, when the difference in the 
dislocation positions in DDD simulations between this and the previous time step is 
no more than 10 _ 5 Atdi S * We see that the relative error in the pair density at different 
stages is no more than 1.5%, while the relative error in pair width is roughly 8 %. In 
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(a) t = 0 


(b) t = 1 


(c) t = 2 



(d) t = 5 




(e) t = 10 (f) t = 26.4 


Fig. 12. Snap shots of the pair density obtained from the DDD and the continuum methods at 
t = 0, 1, 2, 5, 10 and 26.4, where t is measured in unit 27r(l — u)L 2 /(rn g p.b 2 ). 


Fig. 12, snap shots of pair density by using the two methods at t = 0, 1, 2, 5, 10 and 
26.4 are shown. 

We also check the efficiency of the continuum model by keeping all other parame¬ 
ters unchanged while increasing the total number of dislocations N. For this purpose, 
we introduce two quantities T con and T& i s , which denote the wall-clock time it takes 
a simulation to reach the steady state by using the continuum and DDD models, 
respectively. Thus T cori /T&i s becomes a measurement of the computational efficiency 
of using the continuum model against its underlying DDD model. The smaller this 
value is, the higher efficiency the continuum model displays. 

The comparison between the two models for different N is shown in Fig. 13. In 




Fig. 13. (a) T con /T,n s provides a measurement to the computational efficiency exhibited by 
the continuum model compared to its underlying DDD model. The smaller this value is, the more 
efficient the continuum model is. (b) The upscaling errors of the pair density <j>' and the pair width 
defined by Eqs. (7.1) and (7.2), respectively, as the systems attain their steady states with various 
N. 
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Fig. 13(a), T con /T dis is found scaling with N at an exponent of roughly —3.25. When 
the total number of dislocation pairs is increased to 500, the time it takes for the 
continuum model to reach the steady state is roughly 1% of that needed for the DDD 
model. This suggests that the continuum model becomes extremely efficient for a 
large N compared to its DDD counterpart. The greater efficiency displayed by the 
continuum model can be attributed to the fact that an increase in N only brings 
up the computational intensity of performing the DDD simulations, as the governing 
equations (5.22a) and (5.22b) for the continuum model are both independent of N. 

To check the accuracy of the continuum model, we also plot Err^,/ and Err^ given 
by Eqs. (7.1) and (7.2), respectively against N in Fig. 13(b). The coarse-graining 
errors (in the interior region) for both quantities drop with an increasing N. When 
N is 500, the coarse graining error of ( in the interior region measured by Eq. (7.2) 
becomes as good as no more than 1%. This is sensible since the continuum model 
is obtained by taking the asymptotic limit as N —> oo. An increased N effectively 
brings down the truncation errors. 

When the rescaled slip plane gap S is small, the simulation can be speeded 
up using the asymptotic solutions to Eq. (5.22a), rather than numerically solving 
Eq. (5.22a) at each time step. In this scenario, the governing equations at the con¬ 
tinuum level can be asymptotically simplified to 


(7.3a) 


C = S-2S 2 r e ° 



S 3 


and 

(7.3b) 


d(j) 

dt s 


■K 2 S 2 (f>" (f>' + £ • 


C <9r e ° x . 


dx 


S _ d-r e ° xt 
2 dy 


d4>\ d(j) 
dx J dx 


In Table 2, the coarse-graining errors for the pair density and the pair width are 
shown with N = 50 and S = 0.1. The upscaling errors are found well controlled 


t 

5 

10 

20 

50 

100 

200 

250 

300 

Err,/ 

0.0060 

0.0064 

0.0068 

0.0074 

0.0130 

0.0208 

0.0221 

0.0227 

Err ? 

0.0179 

0.0181 

0.0184 

0.0189 

0.0188 

0.0185 

0.0184 

0.0184 


Table 2 

The coarse-graining errors of the pair density distribution and the pair width at various time 
slots. Here S = 0.1, t° = 0.5, dr°/dy = 1, N = 50. Here t is measured in unit 2tt(1 — i/)L 2 /(m g p,b 2 ). 


during the simulations. 

8. Conclusion and further discussion. 

8.1. Conclusion. In this paper, we have studied the collective behaviour of 
a row of dislocation dipoles using matched asymptotic analysis. The discrete-to- 
continuum transition is facilitated by the introduction of two field variables, the dis¬ 
location pair density potential (j) and the dislocation pair width £. The equilibrium 
state at the continuum level is governed by Eqs. (5.11a) and (5.11b), while the dy¬ 
namics at the continuum level is given by Eqs (5.22a) and (5.22b). The following 
conclusions are drawn based on the analysis and the numerical implementation to the 
continuum model. 

Dislocation dipoles are found roughly uniformly distributed in the absence of 
applied stress gradients, and to pile up against a lock when a stress gradient is applied. 
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When the externally applied stress is zero on the primary slip plane y = 0, we 
found three possible equilibrium patterns (as shown in Fig. 4), whose stability depends 
on the value of 4>'S, the ratio of the slip plane gap to the pair center spacing. If <f>'S is 
big (condition (6.4) breaks down), non-localised structures (Equilibrium Type II) are 
the stable configurations. When cj>'S falls below the critical value 0.2465, a localised 
equilibrium structure (Equilibrium Type III) emerges. In this scenario, Equilibrium 
Type II becomes unstable and a natural transition to Equilibrium Type III is observed. 

If the externally applied shear stress Tg Xt is non-negligible, two possible equilib¬ 
rium patterns are found and the one with larger pair width value corresponds to the 
stable configuration as suggested by the shaded region in Fig. 11. 

In the continuum limit, the two field variables introduced evolve on different time 
scales. On the faster scale, the dislocation pairs arrange themselves in local struc¬ 
tures to satisfy the leading-order force balance. On the slower scale, the pair density 
evolves driven by the stress gradient, which is a higher-order effect. Consequently, 
the dipole dynamics, if viewed at the continuum level, can be modelled by an equi¬ 
librium equation for £ given by Eq. (5.22a) and an evolution equation for <f> given by 
Eq. (5.22b). All analytical results have been justified through comparison with the 
underlying DDD simulation results. 

8.2. Implication to the formation of PSBs. The finding of a natural tran¬ 
sition between equilibrium configurations of dislocations in this paper may shed light 
on understanding how localised persistent slip band structures emerge within a non- 
localised channel-vein structure in cyclicly loaded crystals. The analytical results in 
§ 6.1.1 suggest that such a transition takes place, when the slip plane spacings drop 
to a certain value such that the quantity equivalent to (f>'S falls below 0.2465. In a 
cyclicly loaded crystal, it is widely recognised that the gaps between slip planes do 
get narrower as a result of the cross-slip motion of the screw segments in the channels 
shown in Fig. 1(a) (see [14, 22]). Nevertheless, the transition in equilibrium patterns 
due to instability found here may not provide a full explanation to the formation of 
PSBs, because the PSB walls consist more likely of several dislocation pairs rather 
than a single pair as indicated by the Equilibrium Type III. 

8.3. Implication to incorporating SSDs into continuum models of plas¬ 
ticity. The approaches used here to separate physical processes according to their 
associated time scales also provide us some hints towards incorporating statistically 
stored dislocations into continuum models of plasticity consistent with the underlying 
discrete dislocation dynamics. Given t the time scale associated with the continuum 
model (termed as the continuum time scale), it has been shown that the mutual ad¬ 
justment within dislocation pairs characterised by the evolution of £ takes place so fast 
that only its steady (equilibrium) state is observable at the continuum time scale. On 
the other hand, the evolution of the pair density potential 4 > takes place so slowly that 
it appears almost unchanged observed at the continuum time scale. Analogously, a 
well-established continuum model of plasticity is expected to be hierarchic in time. It 
should consist of a set of evolution equations for the geometrically necessary disloca¬ 
tions (GNDs) changing at a normal speed accompanied by another set of quasi-static 
equations describing the SSD structures in equilibrium. 
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